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INFLUENCE OF 2D FINITE EL EMENT MODELING ASSUMPTIONS ON DEBONDING PREDICTION 
FOR COMPOSITE SKIN-STIFFENER SPECIMENS SUBJECTED TO TENSION AND BENDING 

RONALD KRUEGER* AND PIERRE J. MINGUEr 


Abstract. The influence of two-dimensional finite element modeling assumptions on the debonding 
prediction for skin-stiffener specimens was investigated. Geometrically nonlinear finite element analyses using 
two-dimensional plane-stress and plane-strain elements as well as three different generalized plane strain type 
approaches were performed. The computed deflections, skin and flange strains, transverse tensile stresses and 
energy release rates were compared to results obtained from three-dimensional simulations. The study showed that 
for strains and energy release rate computations the generalized plane strain assumptions yielded results closest to 
the full three-dimensional analysis For computed transverse tensile stresses the plane stress assumption gave the 
best agreement. Based on this stud> it is recommended that results from plane stress and plane strain models be used 
as upper and lower bounds. The results from generalized plane strain models fall between the results obtained from 
plane stress and plane strain models. Two-dimensional models may also be used to qualitatively evaluate the stress 
distribution in a ply and the variation of energy release rates and mixed mode ratios with delamination length. For 
more accurate predictions, however, a three-dimensional analysis is required. 

Key words, composite materials, delamination, finite element analysis, fracture mechanics 

Subject classification. Structures and Materials 

1. Background. Many composite components in aerospace structures are made of flat or curved panels with 
co-cured or adhesively bonded frames and stiffeners. Previous investigations of the failure of secondary bonded 
structures focused on loading conditions typically experienced by aircraft crown fuselage panels. Tests were 
conducted with specimens cut from a full-size panel to verify the integrity of the bondline between the skin and the 
flange or frame [1]. A simpler and cheaper specimen configuration that would allow detailed observations of the 
failure mechanism at the skin/flange interface was proposed in reference [2]. The investigations focused on the 
failure mechanisms of a bonded skin/flange coupon configuration subjected to tension, three and four-point 
bending, and combined tension/bending loading [3-6]. 

An analytical methodology was also developed to predict the location and orientation of the first transverse 
matrix crack based on the principal transverse tension stress distribution in the off axis plies nearest the bondline in 
the vicinity of the flange tip [7]. Earlier investigations [2,3] indicated that the matrix cracking occurred when the 
maximum principal transverse tensile stress, ct„, normal to the fiber direction, as calculated by 
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reached the transverse tension strength of the material. Here cr 22 and <y 33 denote the normal stresses and r 23 the shear 
stress in the 2-3 plane perpendicular to the fiber direction. 

A fracture mechanics approach was used to investigate delamination onset once the initial crack had formed. 
The initial crack was modeled as a discrete discontinuity at a location suggested by the microscopic investigation. 
The Virtual Crack Closure Technique (VCCT) was used to calculate mixed mode strain energy release rates G\ (due 
to crack opening), G M , (due to in-plane shear) and Gm (due to in-plane scissoring) [8-10]. Computed total strain 
energy release rates, G T , were compared to the critical value, G c , of the material for the appropriate mixed mode 
ratio (G n /G T ) in order to determine the potential for delamination growth [7]. 

During previous studies [2-4, 6,7], two-dimensional models were used because they are preferred by industry 
due to the fact that modeling time, as well as computational time, remains affordable, especially if many different 
configurations have to be analyzed during the initial design phase. However, it is inherent to any two-dimensional 
finite element model that the geometry, boundary conditions and other properties are constant across the entire 
width. The objective of the current study was to evaluate how the results obtained from two-dimensional finite 
element models of the specimen compared to data obtained from full three-dimensional simulations in order to find 
the limitations of using simple two-dimensional models. Investigating the limits and usefulness of any numerical 
model is in general an important engineering science problem that reaches beyond the current analysis of the 
skin/stringer debond specimen. The present investigation is seen as a contribution to the field complementing earlier 
studies which focussed on the calculation of energy release rates using different types of finite element models [11, 
12 ]. 

2. Problem Description. The current study focused on skin-stiffener debonding resulting from buckling of a 
thin-gage composite fuselage structure as described in reference [7]. In that study, the specimens consisted of a 
bonded skin and flange assembly as shown in Figure 1(a). An IM7/8552 graphite/epoxy system was used for both 
the skin and flange. The skin was made of prepreg tape with a measured average ply thickness of h =0.148 mm and 
had a [45/-45/0/-45/45/90/90/-45/45/0/45/-45] lay-up. The flange was made of a plain-weave fabric with a thickness 
of h =0.212 mm. The flange lay-up was [45/0/45/0/45/0/45/0/45]f, where the subscript “f 1 denotes fabric, “0” 
represents a 0°-90° fabric ply and “45” represents a 0°-90° fabric ply rotated by 45°. The measured bondline 
thickness averaged 0.178mm. Specimens were 25.4-mm wide and 177.8-mm long. The properties of the 
graphite/epoxy material and the adhesive were measured at Boeing and are part of the standard design database for 
the V-22 tilt-rotor aircraft. Typical properties are summarized in Table 1. 

Four quasi-static tension tests with a gage length of 101.6 mm were performed as shown in Figure 1(b) [7]. 
The value of the damage onset load was averaged from four tests and determined to be P=17.8kN with a 
coefficient of variation of 8.9%. Three-point bending tests were performed with a bottom support span of 101 .6 mm 
as shown in Figure 1(c). The value of the damage onset load was averaged from four tests and determined to be 
Q =427.6 kN with a coefficient of variation of 12.8%. The tests were terminated when the flange debonded from 
the skin. Damage was documented from photographs of the polished specimen edges at each of the four flange 
comers identified in Figure 1(a). Comers 1 and 4 and comers 2 and 3 showed identical damage patterns for both 
tests. The damage at comers 2 and 3, formed first and consisted of a matrix crack in the 45° skin surface ply and a 



delamination at the +45°/-45° interface as shown in Figure 2. Therefore, this damage pattern was the focus of the 
current and related earlier analyses [7]. 

The complex nature of the failure observed during the experiments, where the delamination changed across 
the specimen width from a delamination running at the skin surface 45%45° layer interface to a delamination 
propagating in the bondline above, suggests the need for a three-dimensional model [7]. Since many layers of brick 
elements through the thickness would be required to model the individual plies, the size of three-dimensional finite 
element models, however, may become prohibitively large. Two-dimensional models of a longitudinal cut through 
the specimen allow for a detailed modeling of the individual plies and the adhesive in thickness direction. Two- 
dimensional models are preferred by industry due to the fact that modeling time, as well as computational time, 
remains affordable, especially if many different configurations have to be analyzed during the initial design phase. 
These models have been used extensively during previous studies [2-4, 6,7]. However, it is inherent to any two- 
dimensional finite element model that the geometry, boundary conditions and other properties are constant across 
the entire width. 

The current study presents an intermediate step where three-dimensional models were created by extruding 
two-dimensional models across the width. The fact that the delamination changed across the specimen width from a 
delamination running at the skin surface 45°/-45° layer interface to a delamination propagating in the bondline 
above, however, is still not accounted for in this model. Nevertheless, the three-dimensional model takes width 
effects into account and therefore provides insight into the limitations of the use of two-dimensional finite element 
models. For future detailed modeling and analysis of the damage observed during the experiments, the shell/3 D 
modeling technique [13] offers great potential for saving modeling and computational effort because only a 
relatively small section in the vicinity of the delamination front needs to modeled with solid elements. 

The objective of the current study was to evaluate how the assumptions made in developing two-dimensional 
finite element models of the specimen would effect the damage onset prediction for skin-stiffener debonding. 
Geometrically nonlinear finite element analyses were performed using ABAQUS* two-dimensional plane-stress, 
plane-strain and generalized plane strain elements. Additionally, two different methods to create generalized plane 
strain conditions were studied [11,14]. The results were compared to data obtained from full three-dimensional 
simulations to study the feasibility and limitations of using simple two-dimensional models. Computed deflections, 
transverse tensile stresses in the skin surface 45° layer, flange and skin strains, as well as mixed mode strain energy 
release rates were compared. 

3. Analysis Formulation 

3.1. Two-Dimensional Finite Clement Models. Several assumptions are possible in two-dimensional 
models. A plane-stress model imposes the out of plane stresses to be zero (a zz =T xz = r v -=0) at the free edge of each 
ply and allows the displacement to be the free parameter. The plane-strain model, on the other hand, imposes the out 
of plane strains to be zero (£ z: =y xz =y yz = 0) at the free edge, which excessively constrains the plies. Hence, the two 
dimensional plane stress and plane strain conditions may serve as upper and lower limits compared to the three- 
dimensional solution. 


3 



An outline of the two-dimensional model of the specimen showing the boundary conditions, and loads 
applied is shown in Figure 3. Finite element models for an undamaged and a damaged specimen were developed 
using a refined mesh in the critical area of the 45° skin ply where cracking was observed during the tests. Outside 
the mesh refinement area, all plies were modeled with one element through the ply thickness. In the refined region, 
four elements were used per ply thickness for the first two individual flange plies above the bondline and the top 
skin ply below the bondline as shown in Figures 3(b) and (c) and described in more detail in reference [7]. Four 
elements through the thickness were also used to model the adhesive film. Ply properties and adhesive material 
properties used in this study are summarized in Table 1 . 

Damage in the form of a matrix crack in the top 45° skin ply that developed into a delamination in the 45°/- 
45° interface is shown in Figure 2. Delaminations of various lengths were modeled by discrete discontinuities as 
shown in Figure 3(c). The delamination propagating in a bimaterial interface (45°/-45°) may, for small element 
lengths A a at the crack tip, result in incorrect mode separation of the total energy release rate. This effect is caused 
by the oscillatory singularity at the crack tip [15]. To avoid the effect, for the current investigation the element 
length A a was chosen to be about 1/3 of the ply thickness. At the opposite taper, the original mesh created for the 
model of the undamaged specimen was used. Finite element solutions were obtained using the commercial 
ABAQ US “/Standard finite element software. Eight-noded quadrilateral plane-stress (CPS8R) and plane-strain 
(CPE8R) elements with quadratic shape functions and a reduced (2x2) integration scheme were utilized for the 
geometric nonlinear analyses [16]. 

3.2. Generalized Plane Strain Models. One alternative to a full three-dimensional simulation that requires 
about the same modeling effort as the simple two-dimensional models is a generalized plane-strain model. 
Generalized plane strain elements provided by ABAQUS* code are typically used to model a section of a long 
structure that is free to expand or is subjected to loading perpendicular to the plane of modeling. The formulation 
involves a model that lies between two planes that can move with respect to each other, and hence, cause strain in 
the direction perpendicular to the plane of the model that varies linearly with respect to position in the planes [16]. 
Ten-noded quadrilateral generalized plane-strain (CGPE10R) elements with quadratic shape functions and a 
reduced (2x2) integration scheme were used in the current models. The models shown in Figure 2 were modified to 
use generalized plane strain analysis. Each element was assigned two additional nodes, which determined the 
position of the bounding planes and were common to all the elements. 

A different approach to create a generalized plane strain case was suggested by Minguet [14]. The 
generalized plane-strain case is regarded as an intermediate state between plane strain and plane stress, and assumes 
that the out of plane normal strain £-z~Vl£xx* where v L is the laminate Poisson’s ratio, and the shear strain Yxz~Y\z~0. 
With these assumptions, ply stiffnesses were calculated for each ply angle in the specimen and entered as material 
properties for the two-dimensional models shown in Figure 3. Details are given in appendix A. Eight-noded 
quadrilateral plane-stress (CPS8R) elements with quadratic shape functions and a reduced (2x2) integration scheme 
were used for this analysis. 

Another possibility to create a two-dimensional generalized plane strain model was described by Konig [11]. 
The model is obtained by using one row of three-dimensional brick elements instead of plate, shell or membrane 
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elements. Therefore, the two-dimensional models shown in Figures 3(b) and (c) were extruded into the models 
shown in Figures 4(b) and (c). This resulted in a 0.1 mm wide model made of one row of brick elements (C3D20R) 
with quadratic shape functions and a reduced integration scheme. The advantage of this model is the inclusion of 
three-dimensional effects in this two-dimensional finite element model. However, the number of degrees of freedom 
almost quadruple in comparison with membrane elements. The model may be compared to a wall, where one side of 
the model represents the midplane and the other side is a parallel plane. The midplane of the specimen as shown 
Figure 4(b) is a plane of symmetry and remains plane under the applied load. The displacements perpendicular to 
the plane are surpressed (vv=0) at all nodes in the midplane. All displacements in the z-direction at the nodes of the 
parallel plane are coupled to a constant value vr. In ABAQUS H the coupling is achieved by grouping all nodes in the 
parallel plane except for one. The w displacement of all the nodes grouped is then set equal to the vr displacement of 
the one node not included in the group using a linear constraint equation. These coupled displacements at the 
parallel plane will generate stress lesultants at each node of this plane, but the summation of these stress resultants 
has to vanish because of the zero load in this direction. 

3.3. Three-Dimensional Finite Element Model. The three-dimensional model of the specimen with load 
and boundary conditions is shown in Figure 5 for the tension and in Figure 6 for the three-point bending load case. 
The specimen was modeled with ABAQUS K solid twenty-noded hexahedral elements C3D20R with quadratic 
shape functions and a reduced integration scheme. A two-dimensional mesh was made first in the x-y plane, which 
was then extruded across the width to create the model shown in Figure 5. 

The same refined mesh as in the two-dimensional model (Figure 3) was used in the critical area of the 45° 
skin ply where cracking was observed during the tests. Outside the refined area, the mesh was modified to prevent 
the three-dimensional model from becoming excessively large Figure 5(b). The skin plies were grouped into three 
layered elements with -45/45/90, 90/-45/45 and 0/45/-45 respectively, thus taking advantage of the composite solid 
element option in ABAQUS" [16]. The fabric layers were grouped into two layered elements as shown in 
Figure 5(b). In the transition region at the flange several plies were modeled by one element with material properties 
smeared using the rule of mixtures [17]. This procedure used did not ensure the full A-B-D contribution of the plies, 
however, appeared suitable for the small transition region to enforce a reasonable model size. 

For the stress analysis the two-dimensional model was extruded into twenty elements across the width of the 
specimen, with a refined zone (0.74 mm, five elements) near the free edges (z=0.0 mm and z=25.4 mm) as shown in 
Figures 5(a) and (b). This resulted in a model with 64,180 elements and 279,842 nodes yielding at total of 839,526 
degrees of freedom. The nonlinear analysis required about 50 hours of CPU time on a SGI Origin 3200 workstation. 
The fracture mechanics approach requires multiple analyses to calculate the energy release rate for many different 
delamination lengths to obtain a distribution plot. Due to the large computation time required for each analysis it 
was decided to use a smaller model. The two-dimensional model was extruded into only ten elements across the 
width of the specimen as shown in Figures 5(6) and 6. The final model had 32,090 elements, 145,432 nodes 
yielding at total of 436,296 degrees of freedom and required only about eight hours of CPU time on the same 
workstation. This three-dimensional mesh, however, is not fine enough in the vicinity of the free edges (z=0.0 mm 
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and z- 25.4 mm) to accurately capture the influence of the free edges on the distribution of the energy release rates 
across the width. 

For future detailed modeling and analysis of the damage observed during the experiments, the shell/3D 
modeling technique [13] offers great potential for saving modeling and computational effort because only a 
relatively small section in the vicinity of the delamination front needs to modeled with solid elements. However, the 
applicability of the shell/3D modeling technique to the skin/stiffener debond problem needs to be assessed first. The 
scope of the current paper therefore was limited to the comparison of results obtained from different two- 
dimensional models to data obtained from full three-dimensional analysis and resulting recommendations. 

4. Analytical Investigation 

4.1. Global Response. First, the global response of the specimens was computed at the mean quasi-static 
damage onset load determined from experiments. The load-displacement and the load-strain behavior computed 
from different FE models were compared to the corresponding experimental results. This global response was used 
to examine whether the FE models, the boundary conditions, the loads and the material properties used in the model 
yielded reasonable results. Strains were averaged from computed nodal point values over a length corresponding to 
the dimensions of the strain gages shown in Figure 1(a) [7]. 

A schematic of the deformed geometry, the boundary conditions, and the load applied in the simulations is 
shown in Figure 3(a) for the tension load case. In the schematic, the elongation of the specimen caused by the 
applied tensile load is shown. The bending effect caused by the load eccentricity in the flange region, the 
asymmetric layup with respect to the neutral axis, and the membrane stiffening effect is also shown. For the 
specimen subjected to three-point bending the deformed three-dimensional model, the boundary conditions, and the 
load applied in the simulations are shown in Figure 6(a). A detail of the modeled delaminated region is shown in 
Figure 6(b). For short delamination lengths (a < 1.0 mm) crack opening mode 1 was apparent. For longer 
delaminations mode 1 ceased and the surfaces started to overdose as shown in the detail of Figure 6(c). This 
phenomenon was not observed in the results from two-dimensional analyses. 

4.1.1 Tension Test. The load versus displacement plot in Figure 7 shows that plane-strain model exhibited a 
stiffer behavior than the experiments yielding an upper bound while the plane-stress models were more compliant 
yielding a lower bound. The results from two of the generalized plain strain cases and the full three-dimensional 
model fall between the results from the plane stress and plane strain and agree best with the experimental results. 
The generalized plain strain model suggested by Konig studied [11] yields results similar to the plane stress model. 

A comparison of measured strains at the surface of the flange (see Figure 1(a)) and computed results is 
shown in Figure 8. The plane-strain model was stiffer, yielding an upper bound, and the plane-stress model was 
more compliant, yielding a lower bound compared to the results from the other models. The results from the 
generalized plain strain models and the full three-dimensional model agree best with the strain recorded during the 
tests. 

In Figure 9, measured strains at the surface of the top 45° skin ply near the flange tip (see Figure 1(a)) and 
computed surface strains were compared. Scatter in the experiments can not be shown as only one test was equipped 
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with a strain gage on the skin. As mentioned above, the plane-strain models showed a stiffer behavior and the plane- 
stress models a more compliant behavior compared to the tested specimen. The results from two of the generalized 
plain strain cases and the full three-dimensional model fall between the results from the plane stress and plane strain 
simulation, however, appear stiffer when compared to this single test. The generalized plain strain model suggested 
by Konig [11] yields results similar to the plane stress model. 

4.1.2 Three-Point Bending Test. The load versus displacement plot in Figure 10 shows that plane-strain 
model exhibited a stiffer behavior than the experiments (DCDT) yielding an upper bound while the plane-stress 
models were more compliant yielding a lower bound. The results from two of the generalized plain strain cases and 
the full three-dimensional model fall between the results from the plane stress and plane strain and agree best with 
the experimental results. The generalized plain strain model using the ABAQUS* CGPE10R elements yields results 
similar to the plane strain model. 

A comparison of measured slrains at the surface of the flange (see Figure 1(a)) and computed results is 
shown in Figure 1 1. The plane-strain model was stiffer yielding an upper bound, and the plane-stress models more 
compliant yielding a lower bound compared to the results from the other models. The results from the generalized 
plain strain models and the full three-dimensional model fall between these two results, however are consistently 
lower than the experimental results. 

In Figure 12, measured strains at the surface of the top 45° skin ply near the flange tip (see Figure 1(a)) and 
computed surface strains were compared. As mentioned above, the plane-strain model showed a stiffer behavior and 
the plane-stress model a more compliant behavior compared to the tested specimen. The results from two of the 
generalized plain strain cases and the full three-dimensional model fall between the results from the plane stress and 
plane strain simulation. The generalized plain strain model suggested by Konig [11] yields results closer to the plane 
stress model. The generalized plain strain model using the ABAQUS' CGPE10R elements yields results similar to 
the plane strain model. 

In the global response analysis, the load-displacement and the strain-load behavior computed for the tension 
and the three-point bending load cases were compared to the corresponding experimental results. The slightly stiffer 
response of the full three-dimensional model numerical model may be explained by the fact that the material data 
used in the FE simulation originate from the literature. For a consistent simulation, material data should be taken 
from the batch of material that was used to manufacture the specimens. Results from the plane-strain analysis 
indicated a stiffer behavior as expected. This is caused by the constraints inherent to the plane-strain model, 
particularly in the ±45° plies. The plane-stress model which imposes the out of plane stresses to be zero and allows 
the displacement to be the free parameter exhibits - as expected - a more compliant behavior. The stiffnesses of the 
generalized plain strain models fall between the stiffnesses from the plane stress and plane strain assumptions. 

4.2. Local Response 

4.2.1 Stress Analysis. A stress analysis was used to study the initial damage in the form of matrix cracking. 
The maximum principal tensile stress, cr„, normal to the fiber direction was calculated for the static failure load 
using the two-dimensional and three-dimensional finite element models introduced in the previous section. The 


7 



stress distribution in the top 45° skin ply is plotted in the immediate vicinity of the flange tip in Figures 13 and 14 
for the tension and bending tests respectively. In the graph, x=0 mm corresponds to left grip as shown in Figures 2 
to 4, and x=26.4 mm corresponds to the left flange tip. During the fatigue tests, first matrix cracking was observed 
at locations right next to the flange tip which corresponds to the location where peak stresses were calculated [7]. 

The transverse tensile strength and scatter band, for IM7/8552, obtained from three-point bending tests of 
90° lamina, were added to the plot for comparison [18]. Stresses at the static failure load computed from the plane- 
stress analysis correspond with the transverse tensile strength. Stresses from the plane-strain analysis are 
excessively high. This is caused by the constraints inherent to the plane-strain model, particularly in the ±45° plies. 
Plots of the three-dimensional stress distribution across the width of the specimen are shown in Figures 15 and 16. 
The stress distribution has a maximum at jc= 26.6 mm with a peak in the center of the specimen (z=12.7 mm). The 
computed stress drops off towards the edges before it sharply increases near the free edges (z=0.0 mm and z=25.4 
mm). A locally even more refined mesh would be required to capture further details near the edges. 

4.2.2 Fractures Mechanics Analysis. A fracture mechanics approach was used to investigate delamination 
onset once the initial crack had formed. During a series of nonlinear finite element analyses, strain energy release 
rates were computed at each front location for the loads applied in the experiments. A critical energy release rate, 
G c , needs to be determined to predict delamination onset. This critical G is generally identified based on the shape 
of the total energy release rate versus delamination length curve, which is determined through analysis as shown in 
Figures 17 and 18. The G r versus x curve reached a peak at some virtual delamination length and then decreased. 
The delamination was extended to a total simulated length of 2.2 mm to ascertain that the peak value had been 
captured. 

The total energy release rates computed for all two-dimensional and generalized plane strain models are 
plotted in Figures 17 and 18. The values obtained from three-dimensional analysis along the centerline of the 
specimen (z=12.7 mm) was included in Figures 17 and 18. Qualitatively, all results follow the same trend. After a 
small initial drop the computed total energy release rate increases sharply with delamination length, reaches a peak 
value and gradually decreases. As expected the values from plane stress and plane strain analysis form an upper and 
lower bound, except for very short delamination lengths. Results from generalized plane strain models and the 
three-dimensional model fall between the results from two-dimensional analysis. The generalized plain strain model 
suggested by Konig [11] yields results closer to those obtained from the plane stress model. All generalized plane 
strain peak G r results are within 8.5% of the results obtained from full three-dimensional analysis at the center of 
the specimen. 

Unlike the tension loading case, the G r versus x curves are continuously increasing, indicating unstable 
delamination growth. Qualitatively all results follow the same trend. As before the values from plane strain analysis 
form an upper, except for very short delamination lengths. The generalized plain strain model using the ABAQUS® 
CGPE10R elements yields results similar to the plane strain model and yields a lower bound. For larger 
delamination lengths ( a > 0.5 mm) the other generalized plane strain models yield G T results are within 6% of the 
results obtained from full three-dimensional analysis at the center of the specimen. 
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Three-dimensional plots of the distribution of the energy release rate across the width of the specimen are 
shown in Figures 19 and 20. Values at the free edge (z=0.0 mm and z=25.4 mm) have been excluded from the plots 
as the model was not fine enough to accurately capture the influence of the free edges on the distribution of the 
energy release rates. Along the le ngth (x-coordinate) it was observed for the tension case that after a small initial 
drop the computed total energy release rate increases sharply with delamination length, reaches a peak values and 
gradually decreases. Across the width (z-coordinate) the computed total energy release rate gradually increases with 
z before it drops off near the free edge (z=25.4 mm). Unlike the tension loading case, Gj increases continuously 
along the length (x-coordinate), except for the zones near the free edges (z=0.0 mm and z=25.4 mm) where a sharp 
decrease is observed caused by the overclosure of the interfaces after about 1 mm of delamination growth. Across 
the width (z-coordinate) the competed total energy release rate gradually increases with z before it drops off near the 
free edge (z=25.4 mm). 

The variation of mixed mode ratio G s /G r with delamination length is shown in Figures 21 and 22. Here G s 
denotes the sum of the in-plane shearing components G//+G///, and G r denotes the total energy release rate 
G/+G//+G///, where G { is the opening mode. For two-dimensional analyses, where G /n = 0, this definition is equal to 
the previously used definition of the mixed mode ratio, G U /G T . For three-dimensional analysis, which also yields 
results for the scissoring mode G///. the modified definition of G s is introduced. 

For the tension load all analyses yield the same trend, where the delamination initially starts with high 
shearing components, followed by a drop which is equivalent to an increase in G/. For longer delaminations a 
gradual increase in shearing components is observed. For a more realistic comparison of the G }i contribution from 
two-dimensional and three-dimensional models results from the three-dimensional analysis were calculated 
assuming Gs=G// and ignoring G///. As Figure 21 shows, these results fall between the results from two-dimensional 
analysis, for longer delamination lengths. For Gs— G//+G///, the contribution of the in-plane shearing component is 
larger compared to the mode ratios obtained from two-dimensional analysis. 

For the three-point bending load case all analyses yield the same trend as before, where the delamination 
initially starts with high shearing components, followed by a drop which is equivalent to an increase in opening 
mode 1. For longer delaminations the shearing component remains constant. For the case where G^G /f was chosen, 
the results from three-dimensional analysis in the center of the specimen (z=12.7 mm) show a higher shearing 
component compared to the results from two-dimensional analysis. For the other case (Gs=G//+G//) the contribution 
of the shearing component is even larger compared to the mode ratios obtained from two-dimensional analysis. This 
may be an effect caused by the overclosure of the delaminated surfaces observed near the edges in the three- 
dimensional model. A detailed study of this phenomenon would require contact analysis and preferably a refined 
model near the free edges. 

Three-dimensional plots of the distribution of the energy release rate across the width of the specimen are 
shown in Figures 23 and 24. Values at the free edge (z=0.0 mm and z=25.4 mm) have been excluded from the plots 
as the model was not fine enough to accurately capture the influence of the free edges on the distribution of the 
energy release rates. Along the length (x-coordinate) it was observed for the tension case that the delamination 
initially starts with high shearing components, followed by a drop which is equivalent to an increase in opening 
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mode 1. For longer delaminations a gradual increase in shearing components is observed. Across the width (z- 
coordinate) the computed mixed mode ratio drops slightly. For the three-point bending load case it was observed 
along the length (x-coordinate) that the delamination initially starts with high shearing components, followed by a 
drop which is equivalent to an increase in opening mode I. For longer delaminations the shearing component 
remains constant. Unlike the tension loading case, the mixed mode ratio rises sharply at the edges for delaminations 
longer than 1 mm. This is caused by the overclosure of the delaminated surface near the edges where the crack 
opening mode I disappears. 

5. Concluding Remarks. The influence of the assumptions made in developing two-dimensional finite 
element models on skin-stiffener debonding specimens was studied. Geometrically nonlinear finite element analyses 
using two-dimensional plane-stress and plane-strain elements as well as three different generalized plane strain type 
approaches were performed. The computed skin and flange strains, transverse tensile stresses and energy release 
rates were compared to results obtained from three-dimensional simulations. 

The plane stress and plane strain models provided results for skin and flange strains, as well as energy 
release rates, which form an upper and lower bound of the results obtained from full three-dimensional analysis. 
The results from generalized plane strain models fall between the results obtained from plane stress and plane strain 
models. The generalized plane strain models capture the surface strains on the flange and skin very well. Computed 
energy release rates are within 9% of the results from three-dimensional analysis, however, there is not consistency 
across the load cases with respect to which model performs best. Stresses obtained from plane strain and generalized 
plane strain models, were excessively high when compared to results from three-dimensional analysis. The stresses 
from full three-dimensional analysis were lower than results from any of the two-dimensional models, but are 
closest to the plane stress results. 

In general, two-dimensional models are preferred by industry because modeling time, as well as 
computational time, remains affordable, especially if many different configurations have to be analyzed during the 
initial design phase. The effect of two-dimensional modeling assumptions is most marked for 45° plies because of 
their high in-plane Poisson’s ratio, while it is small for 0° and 90°. Based on the results of this investigation - which 
contains 45° tape and fabric plies - it is recommended to use results from plane stress and plane strain models as 
upper and lower bounds. Two-dimensional models may also be used to qualitatively evaluate the stress distribution 
in a ply and the variation of energy release rates and mixed mode ratios with delamination length. The current 
recommendations are based on the analysis of the skin/stringer specimen with the described layup. Further analyses 
are required before the recommendations may be generalized. 
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Appendix A 

Generalized Plane Strain Procedure. The difficulty in using two-dimensional modeling when representing 
laminated composites is that although the laminate may be in a state of plane stress, each lamina is typically not a 
state of plane stress. The effect is most marked for 45° plies because of their high in-plane Poisson’s ratio, while it 
is small for 0° and 90°. There dees not appear to be an exact solution to this problem, however the following 
procedure is an approach that represents a compromise between accuracy and efficiency of two-dimensional 
modeling. 

The traditional three-dimensional stress-strain relationships with the traditional orientations as shown in 
Figure A1 where jc, v, and z are the laminate axes and 1, 2, and 3 the lamina axes can be written as 


<*xx 


£ xx 

%• 


£ yy 

a rr 

£- 

£ zz 

Tyz 


Yyz 

Txz 


Yxz 

Jxv _ 


Yxy . 




a 


where a;,- and i tJ denote the normal and shear stresses and £ ti and y,, the normal and shear strains respectively. In the 
matrix formulation [E] denotes the 6x6 stiffness matrix and [C] the 6x6 compliance matrix. 

The two traditional two-dimensional assumptions are plane strain, where £ vv =y vv -y rr =0 and plane stress, 
where <j xx = r vi .= r v ~=0. An intermediate state is proposed here where is assumed that: £ yv ~- v L e xx and y xy =y yz = 0, where 
v L is a laminate Poisson’s ratio. 

Applying this assumption yields: 


a xx = E \ 1 £ xx ~ v L E I 2 £ xx + E \ 3 £ zz 
a zz = Ei\Sxx ~ Vl E 32 £ xx + E 32 £ zz 
T xz = E 5S Yxz 


or 



and 



where 


H-M" 


The following lamina moduli are then extracted and used as input for a two-dimensional finite element model to be 
run as a plane stress model: 
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The results from the two-dimensional finite element analysis need to be post-processed to recover the 
transverse stresses at the lamina level. 

£yy = —VLE XX 

£-z ~ C2]&xx + Cl2°zz 

a yy “ E 2\ £ xx + E 22 £ yy + ^23 6 zz 

T x v = E 6\ S xx + E 62 £ yy + Ef>2 £ zz 


Once the lamina stresses have recovered, they are rotated into lamina axes. The maximum transverse stress is then 
calculated in the transverse 2-3 plane as shown in Figure A2: 
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TABLE 1. 


Material Properties 


1M7/8552 Unidirectional Graphite/Epoxy Prepreg 

£] i = 161.0 GPa 

£ 22 = H.38 GPa 

£ 33 = 11.38 GPa 

v \2 = 0.32 

v \3 = 0.32 

V 23 = 0.45 

G\2 = 5.17 GPa 

^13 = 5.17 GPa 

G 23 = 3.92 GPa 

IM7/8552 Graphite/Epoxy Plain Weave Fabric 

E if =71.7 GPa 

£22 = 71.7 GPa 

£ 33 = 10.3 GPa 

v\ 2 = 0.04 

vq 3 - 0.35 

V 23 = 0.35 

Gi 2 = 4.48 GPa 

G 13 = 4.14 GPa 

G23= 4.14 GPa 

Grade 5 FM300 Adhesive 

E- 1.72 GPa 

v /= 0.3 

(assumed isotropic) 
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u=v=0 at x=0 

i y 


► X 



v=0 at xsflOl .6 

Static failure load: P=17.8 kN 


undeformed outline 
deformed outline 


(a) Outline of model with load and boundary conditions 



(b) Detail of model of undamaged specimen 



ic ) Detail of deformed model of damaged specimen 
FIGURE 3. Two dimensional finite element model of skin/ flange specimen 
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(c) Detail of deformed model of damaged specimen 
FIGURE 4. Three-dimensional generalized plane strain model of skin! flange specimen 
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(a) Three-dimensional model with load and boundary conditions 



(b) Detail of model of undamaged specimen 



lc ) Detail of deformed model of damaged specimen 
FIGURE 5. Full three-dimensional model of slant flange specimen 


0 at x>101 .6 
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(a) Deformed three-dimensional modi 



(b) Detail of deformed moa 



delamination 

front 


contact 

region 


( c ) Detail of a 

FIGURE 6. Full three-dimensional model of skin/) 
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FIGURE 7. Load-displacement plots for tension tests. 
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FIGURE 13. Computed transverse tensile stress results in top 45° skin ply for tension rest 
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FIGURE 14. Computed transverse tensile stress results in top 45° skin ply for three-point bending test 
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FIGURE 17. Computed total energy release rate for delamination growing from matrix crack 
in skin top 45°/ -45° ply interface for tension test 
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FIGURE 1 8. Computed total energy release rate for delamination growing from matrix crack 
in skin top 45°/-45°ply interface for three-point bending test 




26 




26 


FIGURE 19. Computed total energy release rate distribution for tension test. 



FIGURE 20. C< mputed total energy release rate distribution for three-point bending test. 
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FIGURE 21. Computed mixed mode ratio G JG for delamination growing 
from matrix crack in skin top 45°/-45°ply interface (tension test) 
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FIGURE 22. Computed mixed mode ratio G^/G^for delamination growing 
from matrix crack in skin top 45°/~45°pIy interface (three -point bending test ) 
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FIC URE 23. Computed mixed mode ratio Gs^Gjfor tension test. 
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FIGURE 24. Computed mixed mode ratio C^Gjfor three-point bending test. 
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